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摘要 : 旋转 质心 坐标 系 下 的 椭圆 型 外 限制 性 系 外 行星 三 体 问题 的 哈密 顿 方程 含有 坐标 和 动量 的 
交叉 项 ， 并 且 显 含 时 间 变 量 ， 系 统 不 再 守恒 ， 显 式 力 梯 度 辛 算法 无 法 直接 应 用 。 对 此 ， 通 过 扩大 
相 空 间 将 非 保守 哈密 他 系 统 变换 为 自治 的 哈密 椭 系 统 ， 并 重新 构造 力 梯度 辛 算 法 ， 实 现 力 梯度 辛 
算法 在 椭圆 型 外 限制 性 三 体 问 题 中 的 应 用 。 结 果 表 明 ， 构 造 的 力 梯度 辛 算 法 的 精度 优 于 非 力 梯度 
辛 算法 ， 并 且 优 化 后 的 力 梯 度 辛 算法 的 精度 优 于 未 优化 的 力 梯 度 辛 算法 。 此 外 ， 采 用 优化 的 力 梯 
度 算法 ， 以 及 快速 Lyapunov 指数 对 椭圆 型 外 限制 性 系 外 行星 三 体系 统 进行 相 空 间 扫 描 ， 获 得 了 
各 参数 对 行星 轨道 动力 学 稳定 性 的 影响 。 
关 键 i: 力 梯度 辛 算法 ， 椭圆 型 外 限制 性 三 体 问题 ， 系 外 行星 ;混沌 ， 动 力学 稳定 性 
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三 体 是 指 可 视 为 质点 的 天 体 靠 相互 之 间 万 有 引力 作用 组 成 的 力学 关系 ”， 由 于 一 般 三 体 
问题 是 不 可 积 的 18 维 哈密 顿 系统 ”， 并 且 相 空间 比较 复杂 ， 我 们 常 忽略 第 三 体 的 质量 ， 将 
其 视 为 仅 在 两 主 天 体 的 引力 作用 下 运动 的 试验 粒子 ， 且 不 考虑 它 对 两 主 天 体 的 影响 ， 这 种 
简化 的 三 体 模型 称 为 限制 性 三 体 问题 ”。 由 于 mm « mi,m2， 可 以 不 考虑 行星 m 对 另外 两 主 
天 体 mi 和 m 的 引力 作用 ， 两 主 天 体 m Mm 在 相互 引力 作用 下 做 开 普 勒 运动 ， 行 星 mm 
在 两 个 主 天 体 的 引力 下 运动 。 若 两 主 天 体 m1 和 ms 做 圆 运动 ， 则 为 圆 型 限制 性 三 体 问 题 ; 
若 做 椭圆 运动 ， 则 为 椭圆 型 限制 性 三 体 问 题 。 如 果 行 星 m 的 运行 轨道 在 两 主 天 体 轨道 内 部 ， 
则 为 内 限制 性 三 体 问题 ， 若 在 轨道 外 部 ， 则 为 外 限制 性 三 体 问题 。 
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三 体系 统 中 存在 混沌 运动 ， 混 沌 运动 是 天 体 动力 学 中 一 个 重要 研究 课题 ， 混 沌 意味 着 蝴 
蝶 效 应 ， 即 非 线性 动力 学 系统 对 初始 条 件 的 微小 变化 具有 指数 敏感 性 "“。 了 解 天 体 的 混沌 运 
动 有 助 于 分 析 行 星 的 轨道 演化 和 轨道 稳定 性 ， 而 鉴定 混沌 和 有 序 需 要 可 靠 的 数值 方法 。20 
世纪 80 年 代 初 ， 为 克服 传统 算法 不 足 ，Feng” 与 Ruth” 各 自 独立 提出 能 够 保持 哈密 顿 相 流 
辛 结构 的 数值 积分 算法 ， 称 为 辛 算法 。 这 类 算法 具有 可 使 系统 能 量 误差 无 长 期 变化 的 优点 ， 
通常 称 为 保持 能 量 ， 但 不 是 严格 意义 的 能 量 保持 方法 ”“"。 它 们 特别 适合 用 来 研究 太阳 系 行 
星系 统 、 系 外 行星 系统 和 致密 天 体系 统 的 长 期 动力 学 演化 问题 ”。 辛 算法 可 以 分 为 隐 式 
辛 算法 ”和 显 式 辛 算法 ”。 隐 式 辛 算法 主要 应 用 于 不 可 分 的 哈密 顿 系统 ， 以 隐 式 中 点 法 为 基 
础 构建 高 阶 隐 式 辛 算法 ， 适 用 于 任何 形式 的 哈密 顿 函 数 ， 但 在 积分 过 程 中 由 于 迭代 计算 导致 
计算 效率 低 ， 多 数 情况 下 显 式 辛 算法 通常 要 求 哈密 顿 系统 至 少 分 解 为 两 个 可 积 部 分 。 由 于 显 
式 辛 算法 的 计算 效率 比 隐 式 辛 算法 的 高 ， 当 哈密 顿 系统 能 可 积分 解 时 ， 显 式 辛 算法 应 当 优 先 
考虑 选择 ””。 有 时 需 将 变量 不 可 分 离 哈密 顿 系统 分 解 为 一 个 显 式 求解 部 分 和 另 一 个 隐 式 
求解 部 分 ， 然 后 组 合成 显 隐 混 合 辛 算法 ””， 与 整个 系统 完全 用 隐 式 辛 算法 求解 来 比 也 有 
计算 效率 优势 。 扩 大 相 空间 ， 构 造 显 式 类 辛 算法 “”， 具 有 辛 方法 类 似 性 质 和 计算 效率 快 


TT 


> 的 优点 。 关 于 辛 算法 详细 介绍 可 参阅 雇 新 浩 和 刘 林 ”及 孙 浪 等 人 的 综述 文章 。 

co 1983 4E, Ruth” 针对 具有 动能 了 和 势能 V 两 个 可 积 部 分 的 哈密 顿 系 统 ， 提 出 了 显 式 二 
CO 阶 和 三 阶 辛 算 法 。 随 后 Forest 和 Ruth" 与 Candy 和 Rozmus ™ 的 四 阶 辛 算法 相继 被 提出 ， 
三 1990 Æ, Yoshida 在 Ruth 的 理论 基础 上 用 三 个 二 阶 辛 算法 组 合 得 到 四 阶 显 式 的 Yoshida 
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辛 算法 。1992 Æ, Mclachlan 和 Atela 提出 了 算法 的 优化 概念 ， 将 主要 截断 误差 项 的 影响 
降 到 最 小 ， 得 到 了 优化 的 二 阶 显 式 辛 算法 和 三 阶 Forest-Ruth 辛 算法 。2002 年 ，Omelyan 
等 人 "构造 了 优化 的 四 阶 辛 算法 。 一 个 值得 注意 的 问题 是 Ruth 在 给 出 三 阶 显 式 辛 算法 的 同 
时 还 提出 了 三 阶 力 梯度 算法 ， 就 是 提取 三 阶 截 断 误 差 中 与 力 梯度 有 关 项 加 入 势能 V 来 构建 
算法 ， 亦 即 用 更 改 的 势能 来 构造 辛 算 法 。 后 来 Chin 和 Chen "FE BEE SIE RBS 
四 阶 。Omelyan 等 人 “将 最 优化 思想 引入 到 力 梯 度 辛 算法 中 ， 按 照 对 称 组 合 的 形式 构造 
了 三 阶 和 四 阶 最 优化 力 梯度 辛 算法 ， 并 将 这 两 种 方法 运用 到 天 体力 学 等 领域 中 。 国 内 学 者 也 
对 力 梯 度 辛 算法 的 发 展 和 应 用 拓 广 做 了 一 些 工 作 ”。 实 际 上 ， 力 梯度 辛 算法 的 应 用 不 限 
于 动能 全 为 动量 二 次 标准 形式 。Xu 和 Wu” 发现 了 为 开 普 勒 问题 时 力 梯度 辛 方法 仍 可 应 
í J ERER 将 力 梯度 辛 方法 推广 应 用 于 动能 是 动量 的 二 次 函数 形式 的 圆 型 限制 三 
本 问题 。 近 来 ，Zhang 等 人 ”推广 力 梯度 辛 方法 于 人 了 可 含有 坐标 且 为 动量 二 次 的 形式 。 

研究 发 现 ， 与 同 阶 的 非 力 梯度 辛 算法 相 比 ， 力 梯度 辛 算法 在 开 普 勒 二 体 问 题 、 摄 动 的 开 
普 勒 二 体 问 题 、 分 子 动力 学 多 粒子 问题 、 太 阳 系 N 体 问题 、 量 子 力学 薛 定 翰 方程 以 及 圆 弄 
限制 性 三 体 问 题 模 拟 中 都 比 同 阶 非 力 梯度 辛 算法 好 ， 并 且 优 化 后 的 总 比 未 优化 的 算法 精度 
好 下 加 加 四 网。 旋转 质心 坐标 系 下 的 椭圆 型 外 限制 性 系 外 行星 三 体 问 题 的 哈密 顿 方程 比 
较 复杂 ， 并 且 具 有 坐标 和 动量 的 交叉 项 显 含 时 间 变量 ， 哈 密 顿 量 不 再 守恒 ， 因 此 ， 力 梯度 辛 
算法 的 应 用 遇 到 了 困难 。 然 而 哈密 顿 方程 的 前 半 部 分 具有 解析 解 ， 我 们 通过 扩展 相 空间 ， 将 
时 间 变量 视 为 一 个 新 坐标 ， 可 以 得 到 一 个 新 的 自治 哈密 顿 量 ， 通 过 重新 计算 力 梯度 算 子 的 表 
达 式 进 而 构造 显 式 的 力 梯度 辛 算法 。 本 文 的 主要 目的 是 在 陈云 龙 和 伍 均 号 的 工作 及 Zhang 
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EA 工作 基础 上 进一步 推广 力 梯 度 辛 算法 在 旋转 质心 坐标 系 下 的 椭圆 型 外 限制 性 三 体 问 
题 中 的 应 用 。 在 椭圆 型 外 限制 性 系 外 行星 三 体 问 题 中 采用 通常 的 四 阶 Forest-Ruth 辛 算法 及 
其 优化 算法 ， 与 本 文 构造 的 四 阶 力 梯度 痒 算 法 及 其 优化 算法 进行 数值 模拟 比较 ; 探讨 力 梯度 
辛 算法 在 求解 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 的 精度 和 效率 是 否 优 于 非 力 梯 度 辛 算法 ; 
找到 在 椭圆 型 外 限制 性 系 外 行星 三 体 问 题 中 数值 表现 能 力 最 好 的 积分 器 ， 然 后 应 用 得 到 的 
可 靠 的 混沌 判别 指标 进行 相 空间 扫描 ， 讨 论 在 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 各 参数 
与 行星 轨道 动力 学 稳定 性 的 关系 。 


2 算法 介绍 


假设 动能 T(p) 是 一 个 关于 n 维 动量 p WOK RB BVT = T(p) = p?/2， 势 能 V(g) 是 
n 维 坐标 q 的 函数 。 它 们 组 成 的 哈密 顿 函 数 为 : 


A (p,q) =T(p)+V(q) . (1) 
分 别 以 A 和 B AT AV W Lie BF, WF: 
= 0 
A= T; =— > 2 
3 Ou ( ) 
= 0 
=A ag, (3) 


XE, T; = OT/Op;, V; = 9V/0q;:。 A 算 子 分 别 作 用 于 坐标 和 动量 ， 我 们 可 以 得 到 G = 
Aq, = Tp; I pi = Ap, 三 0， 显 然 4 算 子 仅仅 对 位 置 坐标 有 作用 ， 称 4 为 位 置 型 算 子 。B 
算 子 分 别 作用 于 坐标 和 动量 ， 得 到 Bp; = -Ya 和 Bq = 0，B 算 子 是 动量 型 算 子 。 利 用 两 
个 算 子 对 称 构造 出 二 阶 Verlet 辛 积分 器 : 


eW — exhBehAgghB | (4) 


be T+ V DWN FERFE Eh 是 积分 的 时 间 步 长 ，W 用 Baker-Campbell- 
Hausdroff (BCH) 公式 可 以 写成 : 


44, [B, Al] + 


W =h(A+B) +h! (- | ZIB, [A,B] + o(n°)) (5) 


We HAF [A,B] = AB-BA, F#A C =[B,[A, B]] = [B, AB — BA] = 2BAB- BBA — 
ABB。 容 易 验 证 BBA, = BBAp; =0, ABB, = ABBy, = 0, BAB, =0。 因 此 得 到 三 阶 
RER]: 


af: Š ð 
C = [B, [A, B]] = 2 > ViVi in ee ， (6) 
i i=1 i 


i,j,k=1 
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其 中 ， Vij = OV /Oqi0q;; Vk = OV /04k, Tjk = ƏT /Op; Opp, f = (-Vi, 一 他， —V3, oy , —Va)o 


算 子 C 与 B 属于 同类 动量 型 算 子 ， 称 为 力 梯度 算 子 。 早 在 1983 年 ，Ruth ”在 首次 给 出 三 
阶 显 式 辛 算法 时 ， 将 力 梯 度 算 子 加 入 势能 所 对 应 的 Lie 算 子 中 ， 构 造 了 显 式 的 三 阶 力 梯度 辛 
算法 。 这 种 力 梯度 辛 算法 有 效 地 避免 了 使 用 负 积 分 步 长 ， 相 较 于 同 阶 的 常用 辛 算法 ， 形 式 更 
加 简单 并 且 精 度 更 高 。 后 来 ，Chin 等 人 一 一 ”发 展 了 这 种 力 梯度 辛 算法 并 扩展 到 了 四 
阶 ， 将 其 成 功 应 用 到 经 典 力 学 与 量子 力学 模型 中 ， 证 明了 这 种 力 梯度 辛 算法 的 误差 精度 相 较 
于 同 阶 常用 辛 算法 更 具 优 势 。 显 式 四 阶 力 梯度 辛 算 法 (F4) 为 : 


FA = h(a) MehB HAVEN C MA HB HATS HIYA a 


3 


结合 力 梯度 辛 算法 和 最 优化 辛 算 法 的 思想 ，Omelyan 4AT 在 力 梯度 辛 算 法 中 引入 了 算法 
最 优化 思想 ， 令 五 阶 截 断 误差 各 项 系数 的 平方 和 最 小 ， 得 到 了 最 优化 的 四 阶 力 梯度 辛 算法 
(OF4): 


< 3 = = 3 
OF4= ethB pbha_ (0.5 a)hB+ch Cel 2b)hA (0.5 a)hB+ch C ebhA pah B , (8) 


其 中 算 子 的 系数 为 : 


a = 0.878 936 860 168 070 9 x 10! ， 
b = 0.281 398 061 166 7719 , 
c = 0.306 181 012 236 977 0 x 107? . 


本 文 以 Forest-Ruth 四 阶 辛 算法 (FR)” 及 其 优化 算法 (OFR)™™ 作 参考 ， 评 估 F4 和 
OF4 算法 的 性 能 : 


1 ， (9) 
C=O 


c 1—c 1 一 c c 
wa z hA p(1—2c)hB o3 hA pchB eo ghA 


OFR = "Be 172A hA oXhB pàhA o(1-2(x—8))hB pARA oXhB o 172A hA £hB , (10) 
这 里 各 算 子 的 系数 为 : 


€ = 0.172 086 559 029 5143 , 
A = —0.091 562 030 755 156 78 , 
x = —0.161 621 762 210 7222 . 


正如 引言 所 云 ， 力 梯度 辛 算法 的 应 用 不 限于 动能 了 为 动量 二 次 标准 形式 。 陈 云龙 和 伍 
坎 ” 将 力 梯度 辛 方法 推广 应 用 于 动能 是 动量 的 二 次 函数 形式 的 圆 型 限制 三 体 问题 。 下 面 我 
们 进一步 将 力 梯度 辛 算法 推广 应 用 到 椭圆 型 限制 三 体 问题 中 。 


[Some 
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3.1 ”物理 模型 
在 惯性 坐标 系 中 讨论 椭圆 型 限制 性 三 体 问 题 的 动力 学 积分 ， 可 以 将 哈密 顿 方程 写成 可 
分 离 变量 的 形式 ， 从 而 可 以 利用 显 形式 辛 积 分 器 。 但 由 于 每 一 步 都 要 涉及 两 个 主 天 体位 置 矢 


量 的 计算 ， 而 且 必 须 求 解 开 普 勒 方程 ， 特 别 是 双星 偏心 率 较 大 时 ， 计 算 量 明 显 增加 ， 计 算 效 
率 低 。 而 采用 旋转 质心 坐标 系 时 横 坐 标 始终 在 两 主 天 体 的 连 线 上 ， 两 主 天 体 的 坐标 矢量 是 不 
变 的 ， 可 以 大 大 减少 计算 量 。 因 此 我 们 选择 在 旋转 质心 坐标 系 下 讨论 平面 椭圆 型 限制 性 三 体 


问题 。 
y Al m<m,,m, 
: 、 


记 两 主 天 体 的 质量 为 m, mo， 它们 的 质心 记 作 O. B 
设 以 双星 质心 为 原点 的 质心 坐标 系 为 O 一 En, TBA 
(&,n)， 以 双星 质心 为 原点 的 旋转 质心 坐标 系 为 O 一 xy， 行 星 
坐标 为 (z,y) ( 见 图 1)。zy 面 即 为 两 主 天 体 的 运动 平面 ， 旋 
转 质心 坐标 系 的 旋转 速度 与 两 主 天 体 的 运动 角速度 f 一 致 。 
E ”质心 坐标 系 与 旋转 质心 坐标 系 下 的 位 置 变量 关系 为 : 


m=i-p 
1 ”旋转 质心 坐标 系 O 一 zy 下 T=éc0sf + nsin f 
椭圆 型 限制 三 体 问 题 构 型 y=—ésinf+ncosf 


. (11) 
z=-—ésinf-f+ncosf-f+écosf+7sin f 


y = —€cos f - f — nsin f - f — Ẹsin f + ù cos f 
Ap, f=(i- e)” (1 + e cos f)” 

为 了 简化 问题 ， 取 适当 的 计算 单位 ， 让 引力 常数 G = 1。 质 量 单位 取 m +m = M, 
L= mz/M 为 第 二 天 体 的 无 量 纲 质 量 ， 该 天 体 的 位 置 为 (1 一 ,0)， 男 一 主 天 体 的 无 
ENKEN- u = m/M, MAN (-y,0) KEP MERAH AER r = 
a, (1—e2)/(1+e,cosf) , HP ay, e 和 三 分别 为 ms 相对 于 m 的 轨道 半 长 径 、 轨 道 仿 
心率 和 真 近 点 角 。 

在 上 述 单 位 制 下 ， 广 义 坐 标 g 和 广义 动量 p 分 别 为 : 


q=r 
i. (12) 


qı 一 42 d 
tm r= p kasli n ,9 = gp ANIME RRC: 
0 


[Some 
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1 2 2 
H = 5 (pi + p23) + (gq2p1 — qip2) — 2 (q, f) 
Q (q, f) = : (+8) +U (0) -1ta 
, 1l1+eicosf [2 2 J NAL e 
o a! , (13) 
ry 72 


n= (+u + 


r= (a - ( ) +g 


中 ，el 为 双星 轨道 偏心 率 。 时 含有 位 置 和 动量 的 交叉 项 gop. - apo T 不 是 动 
量 的 严格 二 次 型 ， 并 且 含 有 时 间 变 量 f{， 因 此 ， 哈 密 顿 量 不 守恒 ， 系 统 不 可 分 离 ， 力 梯度 辛 
算法 无 法 直接 使 用 。 
3.2 JARTAK 

我 们 将 f 视 为 华 标 的 一 个 分 量 9o， 并 加 入 对 应 的 动量 pp， 对 式 (13) 进行 相 空 间 扩充 。 
今 po =—-H,w = f， 得 到 一 个 新 的 守恒 的 哈密 顿 量 A* = 0。 哈 密 顿 方程 的 前 半 部 分 是 有 解 
析 解 的 ， 因 此 相 空 间 扩 充 后 的 哈密 顿 方程 仍然 可 以 分 离 成 “动能 T” 和 “势能 VV” 两 个 可 积 
部 分 分 别 求解 。 新 的 哈密 顿 方程 为 


并 


H* = T (q1, q2, P1, P2, Po) + V (q1, gq2, G0) 


1 
T (q1, q2, P1, P2, Po) = 5 (pi+p2) + gp1— qip2+po . (14) 
V (qi, 92,90) = —22 (q1, q2, go) 


动能 了 含有 坐标 和 动量 的 交叉 项 ， 不 再 是 动量 p 的 严格 二 次 型 ， 所 以 4 算 子 的 形式 不 再 是 
A (2) 的 形式 ， 而 是 坐标 和 动量 的 混合 型 算 子 : 


2 ð ð 
A= TR] a 15 
> Pi Og; ne) ( ) 
势能 V 的 Lie 算 子 C 为 : 
C = |B, |A, B]] = [B, AB — BA] =2BAB-— BBA — ABB . (16) 
可 以 验证 
BBA,, = BBA,, =0 
ABB, = ABB,, =0 . (17) 
BAB,, =0 
但 是 ， Ed 
BAB,, = DDI oma (18) 
Pi Lu Luz j j Opi” 


f 
/二 
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(19) 


根据 式 (14) 中 的 T (1; q2, P1, P2, Po) 可 知 ， Xj Æ k AY, Tjk = 0; Hj =k=0 时 ， Tjk = 0; 


显然 ， A J j KHAO, Tjk = 1. 因此 ， 
2 2 9 
0 


ð 
2 (Voi Va + Voz V2) Dns +2 (VaV +V 
0 


上 式 表明 算 子 C 是 两 主 天 体 引 力 的 力 梯度 算 子 。 所 以 ， 旋 转 质心 坐标 系 下 的 李 


性 三 体 问 题 仍然 可 以 用 力 梯 度 辛 算法 求解 。 
接 下 来 ， 我 们 分 别 给 出 该 系统 全 部 分 和 了 
(1) 动能 他 部 分 对 应 的 正则 方程 
利用 式 (14) 可 以 得 到 动能 全 中 算 子 4 (EA 


可 以 得 到 ; 


ð ð 
2V2) =— + 2 (Va Vi + Vo2V2) 二 一 . 


P= pa ~ 
__ or 
d1 Op, = 
or 
d2 Ape 
by coe 
Pı Ba, 
ee 
P2 Bas 


从 积分 初始 状态 量 (goo), q1(0), 42(0)s P1(0)s P2(0)) 出 发 ， 经 过 时 间 /后 获得 该 方程 组 的 解 : 


qo = Go) + f 
qı = U0) cos f + q2(0) Sin f + 


q2 = 42(0) cos f — d1(0) sin f + 
p1 = Pio) cos f + pacoy sin f 
P2 = Pao) COS f- P1(0) sin f 


上 式 意 味 着 混合 算 子 4 是 可 积 可 解 的 。 
(2) 势能 V 部 分 对 应 的 正则 方程 


an an (20) 
圆 型 外 限制 
部 分 对 应 的 正则 方程 。 
日 于 坐标 qi 和 动量 pi 的 微分 方程 组 为 : 
1 
Pı 十 92 
修一 而 (21) 
=p 
= -pı 
f (pro) cos f + p20) sin f) 
f (pac) cos f — pico) sin f) (22) 
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st (14) 可 得 势能 V 中 算 子 B 作用 于 动量 p; 的 微分 方程 组 ， 可 得 : 
i OV 
Po = ~ Ogo = — Vo (q1, 92, qo) 
OV 
ij 二 一 一 = -Vi (q, q2, . 23 
Pı an 1(91 q2 qo) ( ) 
OV 
2 二 一 一 二 一 ) d2, 
p2 Ago 2(qı q2 qo) 
势能 V 中 算 子 C 作用 于 动量 p; 的 微分 方程 组 ， 可 得 : 
Po = 2Vo1(q1, 92, qo) Vi (41; 42, qo) Tr 2V02(q1, 92, qo) V2(q1, q2, q0) 
D1 = 2V1 (q1, q2, qo) Vi (q1, q2, Go) + 2Vi2(q1, 92, qo) V2 (41,9290) > (24) 
P2 = 2V21 (q1, q2, qo) Vi (q1, q2, Go) + 2V22(q1, G2, Go) V2 (q1, q2, Go) 
e 势能 V 部 分 的 解 为 : 
Po = Po(o) — bihVo(q1, G2, qo) + cih? + 2(Vor(G1, G2, qo) Vi (q1, q2, Go) + 
Vo2(q1, q2, qo) V2(q1, 92, do)) 
= pio) — bih Vi (q1, 92, G0) + cih? - 2(Vin (q1, q2, qo) Vi (q1, 92, Go) + (25) 
Vi2(q1, 92, Go) V2 (q1, G2, q0)) 
= D2(0) 一 bihVa(q , 42,90) +h’. 2(V21 (q1, 42, G0) V1 (q , q2, qo) + 


Və2(q1, q2, qo) V2(q1, q2, qo)) 


其 中 ， bi 和 Ci 是 B 算 子 和 C 算 子 各 阶 的 系数 ， (Povo); P10 ), P2(0)) 为 积分 初始 状态 量 


Ha 
[e] 


4 算法 的 性 能 检验 与 应 用 


4.1 ”初始 条 件 

Oo 本 文 以 双星 质量 比 j = 0.1 的 双星 构 型 为 例 ， 研究 椭圆 型 外 限 制 性 系 外 行星 三 体 问题 

双星 的 轨道 半 长 径 取 w = 1 AU。1998 年 Holman 和 Wiegert” 的 研究 发 现 ， 在 双星 外 存在 
个 轨道 临界 半 长 径 ， 在 这 个 临界 半 长 径 以 内 行星 是 不 稳定 的 ， 而 在 这 个 临界 半 长 径 以 外 可 

以 长 期 稳定 存在 。 临 界 半 长 径 经 验 公 式 如 下 : 


= (1.60 E 0.04) T (5.10 主 0.05) Êi t (—2.22 + 0.11) e? + (4.12 am 0.09) HT 
(—4.27 + 0.17) e1u + (—5.09 + 0.11) u? + (4.61 0.36) e212 ， (26) 


其 中 ，ae 是 行星 轨道 临界 半 长 径 ，ei Mu DAEN E mb KMN E ME e (定义 为 
m2/ (mi 十 m2)， 其 中 m 是 指 质量 大 的 恒星 ，rna 是 指 质量 小 的 恒星 )。 


ChinaXiv 合 作 期 刊 


372 天 文学 进展 40 卷 


根据 行星 临界 半 长 径 的 经 验 公 式 (26)， 取 行星 的 初始 半 长 径 为 a = 2.54。， 再 通过 行星 
的 轨道 根 数 a, e, i, 2, w, M (分 别 为 行星 轨道 半 长 径 、 仿 心率、 轨道 倾 角 、 升 交点 经 度 、 近 
点 角 距 、 平 近 点 角 ) 与 位 置 和 速度 之 间 的 关系 求解 出 行星 的 初始 位 置 和 动量 。 由 于 我 们 考虑 
的 是 共 面 情况 ， 所 以 初始 值 取 i= Q = w= M = 0°*。 首 先 根据 牛顿 迭代 法 求解 开 普 勒 方程 


E-—esnE=M , (27) 


得 到 偏 近 点 角 E, Pa AA iE A E 表示 的 位 置 和 动量 的 关系 式 中 : 


r=a(cos—E—e)P+avl—e?sinEQ , (28) 
2 2 
v=— sin EP + VI — 2 cos EQ i (29) 
r r 
其 中 , n= V4/a3 是 平均 角速度 ，r = a(l - ecos E) 是 两 主 天 体 之 间 的 相对 距离 ， 矢量 P 
和 Q 表示 为 : 
ha cos (2 cos w — sin (2 sin w cos i 
E P = | sinQcosw+cosQMsinwcosi |, (30) 
sinw sin t 
— cos (2 sinw — sin 2 cos w cosi 
Q= | -sin Rsinw + cos Rcoswcosi |. (31) 
cos W sin t 
表 1 ”双星 和 行星 的 初始 条 件 便 可 得 到 轨道 演化 的 初始 位 置 和 动量 。 表 1 列举 了 数值 
双星 行星 模拟 的 主要 初始 条 件 。 
a, =1 AU aa = 2.5a¢ 4.2 ”数值 模拟 
4 = 0.1 i= Q=w=M =0° 


令 双星 偏心 率 e = 0.1， 再 根据 式 (26)， 取 行星 轨 
道 半 长 径 ag = 6， 使 用 固定 时 间 步 长 h = 0.08。 选 择 两 
个 不 同 的 行星 轨道 初始 条 件 讨论 对 应 的 轨道 稳定 性 : es = 0.05 和 es = 0.48。 

用 Forest-Ruth 四 阶 辛 算法 (FR)、 优 化 的 四 阶 辛 算法 (OFR)、 四 阶 力 梯 度 辛 算法 (F4) 
和 优化 的 四 阶 力 梯 度 辛 算法 (OF4) 分 别 积 分 两 个 轨道 ， 并 采用 自 适 应 步 长 RKF8(9) 方法 提 
供 更 高 精度 的 参考 解决 方案 。 图 2 给 出 了 各 算法 计算 出 的 两 个 轨道 的 能 量 误差 入 *， 这 里 
AH* = H*(f) — H*(0), H. A*(f) Ml H*(0) ZEN TA f 和 0 的 哈密 顿 量 ， 积 分 时 间 是 10 
(本 文中 积分 时 间 和 积分 步 长 都 是 无 量 纲 单位 )。 显 然 ， 不 管 是 在 轨道 1 还 是 轨道 2 中 ， 对 了 
同类 算法 ， 优 化 的 要 比 未 优化 的 辛 算法 精度 好 ， 并 且 本 文 构造 的 力 梯度 辛 算法 比 同 阶 非 力 梯 
度 的 辛 算法 精度 好 。 我 们 发 现 不 同 轨道 下 相同 算法 的 精度 不 同 ， 为 了 更 清楚 地 了 解 行星 偏心 
率 对 算法 精度 的 影响 ， 图 3 a) 绘制 了 算法 给 出 的 能 量 误差 对 行星 偏心 率 的 依赖 关系 。 结 
发 现在 不 同行 星 偏 心率 时 算法 的 精度 不 同 ， 并 且 在 行星 偏心 率 e > 0.55 时 精度 变化 波动 比 
较 大 。 为 了 进一步 了 解 在 行星 偏心 率 较 大 时 算法 的 精度 比较 ， 我 们 选取 行星 偏心 率 es = 0.6 


el 三 0.1 e2 = 0.05, 0.48 


= 
fe = 

` 
($) 


iF 
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的 轨道 ， 给 出 了 算法 的 误差 精度 图 ， 图 3b) 显示 在 行星 偏心 率 较 大 时 ， 短 时 间 的 积分 情况 


下 依 | 


日 是 力 梯 度 辛 算法 的 精度 优 于 非 力 梯度 辛 算法 ， 但 由 于 行星 偏心 率 较 大 导致 长 时 间 的 


积分 后 算法 都 发 生 了 失真 现象 ”这 说 明 在 行星 偏心 率 较 大 时 算法 不 适用 ， 而 在 偏心 率 小 的 
情况 下 结果 会 更 可 靠 。 总 的 来 说 ， 在 算法 未 发 生 失 真 时 ， 构 造 的 力 梯 度 辛 算法 的 精度 是 优 于 
非 力 梯 度 辛 算法 的 ， 并 且 OF4 算法 的 数值 表现 能 力 最 好 。 此 外 ， 我 们 又 选择 了 双星 质量 比 


人 三 0 


法 的 外 


.2 的 双星 构 型 ， 令 初始 值 e1 = 0.05, as = 6.5， 使 用 固定 时 间 步 长 h = 0.08， 绘 制 各 算 
E 量 误差 与 行星 偏心 率 的 关系 图 4， 得 到 相同 的 结论 。 对 于 大 偏心 率 轨道 需 引 入 时 间 变 


换 实施 辛 算法 ”™ ， 这 一 问题 在 未 来 的 工作 中 讨论 。 


lg|AH*| 


lg f 
a) 


YE: 使 用 的 算法 分 别 为 FR, OFR, F4, OF4 和 RKF8(9). e1 = 0.1, aa = 6， 时 间 步 长 是 hh = 0.08， 积 分 时 


间 为 105。 
图 2 ”在 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 不 同 算法 分 别 解决 两 个 轨道 的 能 量 误差 值 


"FR 
。OFR 
^F4 

v OF4 

+ RKF8(9) 


YE: el = 0.1, a2 = 6， 同 一 时 间 步 长 h = 0.08, ASE e2 值 ， 积 分 10° 之 后 得 到 的 误差 。 
a) 在 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 各 算法 的 能 量 误 差 与 行星 偏心 率 的 关系 ; b) 行星 偏心 率 
eo = 0.6 时 各 算法 的 能 量 误差 值 
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YE: 初始 值 e1 = 0.1, as = 6.5， 时 间 步 长 是 hh = 0.08， 积 分 时 间 为 105。 


4 WERE u = 0.2 的 双星 构 型 中 各 算法 的 能 量 误差 与 行星 偏心 率 的 关系 


4.3 ”动力 学 性 质 

由 于 OF4 方法 在 能 量 守恒 和 算法 精度 上 具有 最 好 的 数值 表现 能 力 ，FLI 在 效率 和 灵敏 
度 上 都 优 于 最 大 Lyapunov 指数 ， 因 此 我 们 用 OF4 算法 和 FLI 混沌 指标 来 扫描 由 椭圆 型 外 
限制 性 系 外 行星 三 体系 统 中 的 行星 偏心 率 e。 和 行星 轨道 半 长 径 ay 所 决定 的 参数 空间 ， 研 究 
它们 与 行星 轨道 动力 学 稳定 性 的 关系 。 

确定 轨道 混沌 性 质 的 方法 主要 有 : Poincore 截面 、 最 大 Lyapunov 指数 、 快 速 Lyapunov 
指标 、 较 小 排列 指标 、0-1 指标 和 频谱 分 析 等 ”。Lyapunov 指数 和 快速 Lyapunov 指数 是 
区 分 系统 有 序 混沌 最 常用 的 方法 ， 它 们 克服 了 Poincore 截面 仅 限于 相 空 间 为 四 维 的 保守 系 
统 的 局 限 。Lyapunov 指数 是 衡量 两 邻近 轨道 随时 间 平 均 指 数 分 离 比 的 指标 ， 能 够 反映 混沌 
的 强度 。 通 常 采用 最 大 Lyapunov 指数 ， 其 计算 方法 主要 有 变 分 法 和 两 粒子 法 ”。 构 造 变 分 
方程 及 其 解 通常 不 容易 ， 考 虑 到 应 用 的 方便 ， 一 般 采 用 两 邻近 轨道 的 差 来 代替 变 分 方程 的 
解 ， 即 两 粒子 法 ， 定 义 为 : 


_ 1, d(f) 

oe ay 

HF, d(f) 和 a(0) 是 两 邻近 轨道 分 别 在 f 时 刻 和 0 时 刻 的 距离 。 需 要 注意 的 是 ， 初 始 距 

离 d (0) 的 选择 必须 适当 ， 如 果 dq (0) 太 大 ， 则 相 邻 两 轨道 间 的 差 值 将 与 变 分 方程 的 解 差 别 太 

大 ; 如 果 d (0) 太 小 ， 则 舍 入 误差 会 增 大 。Tancredi 等 人 四 认为 ， 就 双 精 度 而 言 初始 距离 的 

最 佳 选 择 应 为 4(0) = 1078. sh (32) 随时 间 指 数 增 大 ， 会 出 现 溢 出 现象 ， 需 要 重 整 化 处 理 ， 

并 且 需 要 适当 考虑 重 整 化 次 数 。 混 沌 有 序 的 判断 标准 是 : 若 lg| 和 | 趋向 于 一 个 固定 常数 ， 有 
界 系 统 是 混沌 的 ; 若 lg | 和 | 趋向 负 无 穷 ， 则 系统 是 有 序 的 。 

与 最 大 Lyapunov 指数 相 比 ， 快 速 Lyapunov 指数 (FLI) 计算 过 程 更 简单 ， 能 够 更 快 、 

更 灵敏 地 区 分 有 序 轨道 和 混沌 轨道 。 由 于 具有 快速 识别 混沌 的 优点 ，FLI 常 被 用 来 追踪 行 
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(33) 


3 期 REG, SE: 力 梯度 辛 算法 在 外 限制 性 系 外 行星 三 体 问题 中 的 应 用 
星 的 轨道 稳定 性 和 研究 其 他 动力 学 问题 FPL 也 可 以 用 两 粒子 法 计算 ”， 定 义 为 : 
_, d(f) 
FLI = le Ti . 


快速 Lyapunov 指数 也 需要 不 断 重 整 化 ， 并 且 重 整 化 次 数 相对 于 最 大 Lyapunov 指数 少 


=. Saf >l 
相距 d (0) 处 ， 然 后 
时 轨道 有 序 ， 呈 指数 增长 时 
5 和 图 6 分 别 给 出 了 各 算法 求解 两 个 轨道 的 最 大 Lyapunov 指数 和 FLI 示意 图 。 


Ig\A| 


lg 


图 5 


图 


注 


时 实施 重 整 化 ， 即 沿 着 向 量 
再 从 该 点 开始 进行 下 次 积分 。 
了 界 轨道 混沌 。 


主 ， 积 分 时 间 为 105。 


椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 不 同 算法 分 别 求解 两 个 轨道 的 最 大 Lyapunov 指数 


lglM 
I 1 
N = 


(A) 的 方向 将 邻近 轨道 强行 拉 回 到 与 轨道 
混沌 有 序 的 判断 标准 是 : 


5a) 一 e) 中 ， 轨 道 1 的 Lyapunov 指数 随时 间 线 怕 


个 稳定 值 ， 


FE 减少 ， 但 
这 意味 着 轨道 1 是 有 序 的 ;， 男 一 方面 ， 轨 道 2 的 最 大 Lyapunov 指数 在 积分 期 间 


FLI 呈 线 性 增长 


在 积分 期 间 并 不 倾向 于 一 


Ns 
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趋向 一 个 稳定 值 ， 毫 无 疑问 轨道 2 是 混沌 的 。 在 图 6a)—e) 中 ， 轨 道 1 的 FLI SNS 
数 增长 ， 轨 道 是 有 序 的 ; 轨道 2 的 PL 呈 指 数 增长 ， 是 混沌 的 。 可 以 看 到 最 大 Lyapunov 
指数 和 FLI 在 不 同 算法 积分 下 给 出 了 一 致 的 有 序 、 混 沌 判断 结果 : 轨道 1 是 有 序 的 ， 轨 道 
2 是 混沌 的 。 但 是 ， 在 区 分 轨道 混沌 有 序 时 快速 Lyapunov 指数 比 最 大 Lyapunov 指数 更 敏 
感 。 利 用 最 大 Lyapunov 指数 去 区 分 混沌 有 序 至 少 需要 积分 105， 而 ELI 只 需要 积分 104， 
所 以 FLI 是 判断 混沌 更 好 的 指标 选择 。 


FLI 


FLI 


主 ;， 积 分 时 间 为 105。 
图 6 在 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 不 同 算法 分 别 求解 两 个 轨道 的 FLI 示意 图 


ee 


除了 上 面 的 两 条 轨道 ， 我 们 还 模拟 若干 不 同 轨道 下 的 FLI 示意 图 ， 结 果 和 图 6 一 样 ， 
这 里 就 不 再 展示 。 模 拟 结果 表明 ，FLI = 5 可 以 作为 轨道 有 序 或 有 界 混沌 的 临界 值 ， 即 当 
轨道 的 FLI < 5 时 ， 轨 道 是 有 序 的 ， 当 FL > 5 时 轨道 是 有 界 混沌 或 者 无 界 不 稳定 的 。 表 
2 给 出 了 在 图 2a)、 图 5 和 图 6 中 各 算法 求解 轨道 1 所 需 的 CPU 时 间 ， 对 比 了 计算 成 本 。 
毫 无 疑问 ， 力 梯度 辛 算 法 的 计算 需要 一 些 额外 的 计算 成 本 ， 但 不 是 很 昂贵 。 详 细 的 CPU 时 
间 见 表 2。 
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表 2 不同 算法 求解 轨道 1 所 需 的 CPU 时 间 


T CPU 时 间 
Ne Err 入 FLI 
FR 7 90 
OFR 7 102 8 
F4 7 105 8 
OF4 9 126 8 
RKF8(9) 31 496 64 


注 : Err, 入 和 FLI 对 应 的 时 间 分 别 为 求解 能 量 误差 (积分 时 间 为 105)、 最 大 
Lyapunov 指数 (积分 时 间 为 105) 和 FLI (积分 时 间 为 105) 的 CPU 时 间 。 


令 初 始 值 el = 0.1, aa = 6 固定 不 变 ， 将 es 的 初始 值 由 0 以 0.01 的 间隔 递增 到 1, 
每 条 轨道 使 用 OF4 FIERA 104 后 计算 FL 值 ， 同 时 使 用 RKF8(9) 方法 作 参 考 。 图 
7a) 显示 ， 当 行星 偏心 率 ex € (0.0,0.44) 时 ， 轨 道 是 稳定 有 序 的 ， 而 对 于 更 大 的 偏心 率 
ez € (0.45,1) 时 ， 轨 道 是 E ET 稳定 的 。 这 一 结果 表明 ， 在 椭圆 型 外 限制 性 系 
外 行星 三 体系 统 中 ,行星 轨道 动力 学 稳定 性 与 行星 轨道 偏心 率 之 间 存 在 明确 的 依赖 关系 ， 行 
星 偏 心率 越 大 ， 轨 道 就 越 容易 出 现 有 界 混沌 或 者 无 界 不 稳定 。 


到 


ELI 


注 : 其 中 圆 形 的 是 OF4 计算 的 结果 ， 三 角形 的 是 RKF8(9) 计算 的 结果 。 相 同 的 时 间 步 长 h = 0.08, $ 
FLI 是 在 积分 时 间 104 之 后 得 到 的 。 红 色 虚 线 对 应 FLI = 5， 这 是 混沌 区 域 和 有 序 区 域 之 间 的 边界 。 
FLI <5 且 积分 后 的 轨道 半径 R < 200 AU 时 ， 是 有 序 轨道 ; 4 FLI > 5 时 ，R < 200 AU 且 积分 后 的 
心率 ez < 1 时 ， 是 有 界 混沌 轨道 ; 4 R>200 AU 或 es > 工时， 是 无 界 不 稳定 轨道 。 绿 色 的 圆 形 和 
表示 有 序 轨道 ， 黑 色 的 圆 形 和 三 角形 表示 有 界 混沌 轨道 ， 红 色 的 圆 形 和 三 角形 表示 无 界 不 稳定 轨道 。 


7 ”使 用 OF4 和 RKF8(9) 方法 计算 出 的 PLI 分 别 跟 e。 和 a2 的 变化 关系 


S ff lk > 


Seq] 


初始 值 el = 0.1, e2 = 0.05 固定 不 变 ， 将 aa 的 初始 值 由 1 以 0.05 的 间隔 递增 到 7， 对 
不 同 的 行星 轨道 半 长 径 丝 积 jd 104 之 后 计算 FLI 值 ， 从 而 得 到 不 同 的 行星 轨道 半 长 径 as 对 
应 的 轨道 FLI 值 的 变化 ， 得 出 不 同 轨道 的 混沌 有 序 性 。 图 7b) 将 所 有 结果 展现 出 来 ， 可 见 
在 a = 2.3 ita eee 7 od a, 在 这 个 值 以 外 是 有 序 的 ， 而 在 这 个 值 以 


ao 


ChinaXiv 合 作 期 刊 


202306.00387v1 


chinaXiv 


ChinaXiv 合 作 期 刊 


378 天 文学 进展 40 卷 


内 是 混沌 或 者 不 稳定 的 ， 即 az € (2.4,7) 时 ， 轨 道 是 有 序 的 ， 而 在 a € (1,2.3) 时 ， 轨 道 是 
有 界 混沌 或 者 无 界 不 稳定 的 。 显 然 ， 椭 圆 型 外 限制 性 系 外 行星 三 体系 统 的 动力 学 稳定 性 与 行 
星 的 轨道 半 长 径 也 有 直接 联系 ， 行 星 的 轨道 半 长 径 越 小 ， 行 星 轨道 越 容易 发 生 混沌 甚至 不 稳 
定 现 象 。 

为 了 更 好 地 了 解 行星 偏心 率 和 行星 轨道 半 长 径 对 轨道 混沌 边界 的 影响 ， 我 们 令 el = 0.1 
不 变 ， 扫 描 一 个 由 参数 es 和 as 决定 的 二 维 相 空间 ， 探 究 二 维 平面 生成 FLI 的 轮廓 图 。 由 
8 可 见 ， 对 于 椭圆 型 外 限制 性 系 外 行星 三 体系 统 来 说 ， 在 行星 偏心 率 es 偏 大 而 行星 轨道 
半 长 径 as 偏 小 的 区 域 更 容易 出 现 混 沌 甚至 轨道 不 稳定 。 


7 


区 cn D 
FLI 


N 


0.0 


VE: 初始 值 el = 0.1 是 固定 的 。 当 FLI < 5 且 积分 后 的 轨道 半径 尺 < 200 AU 时 ， 是 有 序 轨道 ; 当 
FLI>5, R < 200 AU 且 积 分 后 的 偏心 率 es < 1 时 ， 是 有 界 混沌 轨道 ; 4 R > 200 AU Res S>Im, Æ 
无 界 不 稳定 轨道 。 蓝 色 区 域 表 示 有 序 轨道 ， 绿 色 区 域 表 示 有 界 混沌 轨道 ， 橙 色 区 域 表 示 无 界 不 稳定 轨道 。 


8 ”利用 FLI 扫描 参数 ex 和 a2 的 二 维 空间 来 寻找 混沌 


5 总 结 与 展望 


本 文通 过 相 空间 扩 展 将 非 保守 的 哈密 顿 系统 变换 为 自治 的 哈密 顿 系统 ， 重 新 计算 力 梯 
度 算 子 的 表达 式 并 构造 力 梯 度 辛 算法 ， 在 平面 椭圆 型 外 限制 性 系 外 行星 三 体 问 题 中 实现 了 
力 梯度 辛 算法 的 应 用 。 选 取 不 同 的 轨道 讨论 了 四 种 辛 算法 的 精度 和 稳定 性 。 数 值 模拟 结果 表 
明 ， 在 平面 椭圆 型 外 限制 性 系 外 行星 三 体 问题 中 ， 本 文 构造 的 力 梯 度 辛 算法 的 精度 优 于 非 力 
梯度 辛 算法 ， 并 且 优 化 后 的 辛 算法 的 精度 优 于 未 优化 的 。 优 化 的 四 阶 力 梯 度 辛 算法 具有 最 
高 的 数值 精度 ， 四 种 辛 算法 的 精度 比较 为 OF4 > F4 > OFR > FR。 此 外 ， 本 文采 用 OF4 
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算法 和 快速 Lyapunov 指数 进行 相 空 间 扫 描 ， 研 究 了 椭圆 型 外 限制 性 三 体系 统 中 行星 的 轨道 
偏心 率 和 轨道 半 长 径 变 化 对 行星 轨道 稳定 1 


in 


行星 


行星 越 容易 出 现 混沌 


易 出 
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生 的 影响 。 结 果 表 明 ， 在 椭圆 型 外 限制 性 系 外 行 


现 混 沌 


我 们 在 


星 三 体系 统 中 ， 行 星 轨 道 动力 学 稳定 
的 轨道 1 


T 


性 与 行星 的 轨道 偏心 率 及 轨道 半 长 径 都 有 着 密切 联系 。 
心率 越 大 ， 行 星 越 容易 发 生 混沌 甚至 轨 i 


L 道 不 稳定 现象 ， 行 星 轨道 半 长 径 越 小 ， 


情况 ， 忽 略 了 很 多 影响 


型 外 限制 


t 至 不 稳定 ， 在 行星 偏心 率 人 
甚至 不 稳定 。 
这 里 讨论 的 椭圆 


ij 大 而 行星 轨道 半 长 径 小 的 区 域 ， 轨 道 更 容 


性 三 体 问题 其 实 都 是 假设 了 理想 状态 ， 并 且 讨 论 的 都 
是 因素 。 未 来 希望 对 天 体 运 动 的 研究 不 仅 停留 在 这 种 理想 情况 
下 ， 更 希望 找到 数值 表现 能 力 更 好 的 时 间 变 换 力 梯度 辛 算法 和 混沌 指标 ， 以 解决 空间 椭圆 型 
出 性 三 体 问 题 和 椭圆 型 内 限制 性 三 体 问题 ， 并 进一步 考虑 影响 行星 轨道 运行 的 其 他 因 


素 ， 更 加 准确 地 刻画 天 体 动力 学 的 定性 性 质 。 
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Application of Force Gradient Symplectic Algorithm to 
External-restricted Three-body Problem with two Stars and 


an Extroplanet 


WANG Ya-rut?, LIU Fu-yao!, WANG Ying!?, SUN Wei!?, 
ZHENG Jing-jing'?, XIAO Qian-qian!? 
(1. Shanghai University of Engineering Science, School of Mathematics, Physics and Statistics, Shanghai 


201620, China ; 2. Shanghai University of Engineering Science Center of Applications and Research of 
Computational Physics, Shanghai 201620, China) 


Abstract: The three-body system is not integrable, and can be chaotic in some circum- 
stances. The description of chaotic motion depends on reliable numerical methods and chaos 
indicators. Symplectic algorithms can keep Hamiltonian phase flow and overcome the prob- 
lem that traditional algorithms have secular drifts in energy errors. Because of this, they 
are the best integrators for studying the long-term dynamic evolution of the Hamiltonian 
system. The Hamiltonian of an elliptic external-restricted three-body problem with two s- 
tars and an extroplanet in the rotating centroid coordinate system contains the cross term 
of coordinates and momentum, and explicitly depends on time variable. The Hamiltonian 
system is no longer conserved, and the explicit force gradient symplectic algorithm cannot be 
directly applied. Extending the phase space to transform the non-conservative Hamiltonian 
system into a conservative Hamiltonian system, we can establish force gradient symplectic 
algorithms. The standard fourth-order symplectic algorithm, the optimized fourth-order 
symplectic algorithm, the optimized fourth-order force gradient symplectic algorithm and 
the optimized fourth-order force gradient symplectic algorithm are independently used to 
solve the elliptic external-restricted three-body problem. It is found that the accuracy of 
the force gradient symplectic algorithm constructed in this paper is better than that of the 
non-force gradient symplectic algorithm, and the accuracy of the optimized force gradient 
symplectic algorithm is better than that of the symplectic algorithms without optimization. 
The fourth-order optimized force gradient algorithm and fast Lyapunov index are used to 
scan the phase space of the elliptic external-restricted three-body system, and the effects of 


each parameter on the planetary orbit dynamic stability are obtained. 


Key words: force gradient symplectic algorithm; elliptic external-restricted three-body 


problem; extroplanet; chaos; dynamical stability 


